// This function is ported from libslater library
// https://github.com/nubakery/libslater

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "rys_roots.h"
#include "roots_xw.dat"

/*
from mpmath import *
mp.dps = 25
for i in range(14):
    for j in range(14):
        print ' %.19e' % cos(pi*j*(2*i+1)/28)
*/
static const double COS_14_14[] = {
 1.                       ,
 9.9371220989324260398e-01,
 9.7492791218182361934e-01,
 9.4388333030836757409e-01,
 9.0096886790241914600e-01,
 8.4672419922828412453e-01,
 7.8183148246802980363e-01,
 7.0710678118654757274e-01,
 6.2348980185873348336e-01,
 5.3203207651533657163e-01,
 4.3388373911755812040e-01,
 3.3027906195516709698e-01,
 2.2252093395631439288e-01,
 1.1196447610330785560e-01,
 1.                       ,
 9.4388333030836757409e-01,
 7.8183148246802980363e-01,
 5.3203207651533657163e-01,
 2.2252093395631439288e-01,
-1.1196447610330785560e-01,
-4.3388373911755812040e-01,
-7.0710678118654757274e-01,
-9.0096886790241914600e-01,
-9.9371220989324260398e-01,
-9.7492791218182361934e-01,
-8.4672419922828412453e-01,
-6.2348980185873348336e-01,
-3.3027906195516709698e-01,
 1.                       ,
 8.4672419922828412453e-01,
 4.3388373911755812040e-01,
-1.1196447610330785560e-01,
-6.2348980185873348336e-01,
-9.4388333030836757409e-01,
-9.7492791218182361934e-01,
-7.0710678118654757274e-01,
-2.2252093395631439288e-01,
 3.3027906195516709698e-01,
 7.8183148246802980363e-01,
 9.9371220989324260398e-01,
 9.0096886790241914600e-01,
 5.3203207651533657163e-01,
 1.                       ,
 7.0710678118654757274e-01,
 0.                       ,
-7.0710678118654757274e-01,
-1.                       ,
-7.0710678118654757274e-01,
 0.                       ,
 7.0710678118654757274e-01,
 1.                       ,
 7.0710678118654757274e-01,
 0.                       ,
-7.0710678118654757274e-01,
-1.                       ,
-7.0710678118654757274e-01,
 1.                       ,
 5.3203207651533657163e-01,
-4.3388373911755812040e-01,
-9.9371220989324260398e-01,
-6.2348980185873348336e-01,
 3.3027906195516709698e-01,
 9.7492791218182361934e-01,
 7.0710678118654757274e-01,
-2.2252093395631439288e-01,
-9.4388333030836757409e-01,
-7.8183148246802980363e-01,
 1.1196447610330785560e-01,
 9.0096886790241914600e-01,
 8.4672419922828412453e-01,
 1.                       ,
 3.3027906195516709698e-01,
-7.8183148246802980363e-01,
-8.4672419922828412453e-01,
 2.2252093395631439288e-01,
 9.9371220989324260398e-01,
 4.3388373911755812040e-01,
-7.0710678118654757274e-01,
-9.0096886790241914600e-01,
 1.1196447610330785560e-01,
 9.7492791218182361934e-01,
 5.3203207651533657163e-01,
-6.2348980185873348336e-01,
-9.4388333030836757409e-01,
 1.                       ,
 1.1196447610330785560e-01,
-9.7492791218182361934e-01,
-3.3027906195516709698e-01,
 9.0096886790241914600e-01,
 5.3203207651533657163e-01,
-7.8183148246802980363e-01,
-7.0710678118654757274e-01,
 6.2348980185873348336e-01,
 8.4672419922828412453e-01,
-4.3388373911755812040e-01,
-9.4388333030836757409e-01,
 2.2252093395631439288e-01,
 9.9371220989324260398e-01,
 1.                       ,
-1.1196447610330785560e-01,
-9.7492791218182361934e-01,
 3.3027906195516709698e-01,
 9.0096886790241914600e-01,
-5.3203207651533657163e-01,
-7.8183148246802980363e-01,
 7.0710678118654757274e-01,
 6.2348980185873348336e-01,
-8.4672419922828412453e-01,
-4.3388373911755812040e-01,
 9.4388333030836757409e-01,
 2.2252093395631439288e-01,
-9.9371220989324260398e-01,
 1.                       ,
-3.3027906195516709698e-01,
-7.8183148246802980363e-01,
 8.4672419922828412453e-01,
 2.2252093395631439288e-01,
-9.9371220989324260398e-01,
 4.3388373911755812040e-01,
 7.0710678118654757274e-01,
-9.0096886790241914600e-01,
-1.1196447610330785560e-01,
 9.7492791218182361934e-01,
-5.3203207651533657163e-01,
-6.2348980185873348336e-01,
 9.4388333030836757409e-01,
 1.                       ,
-5.3203207651533657163e-01,
-4.3388373911755812040e-01,
 9.9371220989324260398e-01,
-6.2348980185873348336e-01,
-3.3027906195516709698e-01,
 9.7492791218182361934e-01,
-7.0710678118654757274e-01,
-2.2252093395631439288e-01,
 9.4388333030836757409e-01,
-7.8183148246802980363e-01,
-1.1196447610330785560e-01,
 9.0096886790241914600e-01,
-8.4672419922828412453e-01,
 1.                       ,
-7.0710678118654757274e-01,
 0.                       ,
 7.0710678118654757274e-01,
-1.                       ,
 7.0710678118654757274e-01,
 0.                       ,
-7.0710678118654757274e-01,
 1.                       ,
-7.0710678118654757274e-01,
 0.                       ,
 7.0710678118654757274e-01,
-1.                       ,
 7.0710678118654757274e-01,
 1.                       ,
-8.4672419922828412453e-01,
 4.3388373911755812040e-01,
 1.1196447610330785560e-01,
-6.2348980185873348336e-01,
 9.4388333030836757409e-01,
-9.7492791218182361934e-01,
 7.0710678118654757274e-01,
-2.2252093395631439288e-01,
-3.3027906195516709698e-01,
 7.8183148246802980363e-01,
-9.9371220989324260398e-01,
 9.0096886790241914600e-01,
-5.3203207651533657163e-01,
 1.                       ,
-9.4388333030836757409e-01,
 7.8183148246802980363e-01,
-5.3203207651533657163e-01,
 2.2252093395631439288e-01,
 1.1196447610330785560e-01,
-4.3388373911755812040e-01,
 7.0710678118654757274e-01,
-9.0096886790241914600e-01,
 9.9371220989324260398e-01,
-9.7492791218182361934e-01,
 8.4672419922828412453e-01,
-6.2348980185873348336e-01,
 3.3027906195516709698e-01,
 1.                       ,
-9.9371220989324260398e-01,
 9.7492791218182361934e-01,
-9.4388333030836757409e-01,
 9.0096886790241914600e-01,
-8.4672419922828412453e-01,
 7.8183148246802980363e-01,
-7.0710678118654757274e-01,
 6.2348980185873348336e-01,
-5.3203207651533657163e-01,
 4.3388373911755812040e-01,
-3.3027906195516709698e-01,
 2.2252093395631439288e-01,
-1.1196447610330785560e-01,
};

static void _clenshaw_dc(double *rr, const double *x, double u, FINT nroot)
{
    double d[14];
    double g[14];
    double u2 = u * 2.;
    FINT i, k;

    for (i = 0; i < nroot; ++i) {
        d[0 ] = 0; g[0 ] = x[13+14*0 ];
        d[1 ] = 0; g[1 ] = x[13+14*1 ];
        d[2 ] = 0; g[2 ] = x[13+14*2 ];
        d[3 ] = 0; g[3 ] = x[13+14*3 ];
        for (k = 11; k >= 1; k-=2) {
            d[0 ] = u2 * g[0 ] - d[0 ] + x[k+1+0 *14];
            d[1 ] = u2 * g[1 ] - d[1 ] + x[k+1+1 *14];
            d[2 ] = u2 * g[2 ] - d[2 ] + x[k+1+2 *14];
            d[3 ] = u2 * g[3 ] - d[3 ] + x[k+1+3 *14];
            g[0 ] = u2 * d[0 ] - g[0 ] + x[k+  0 *14];
            g[1 ] = u2 * d[1 ] - g[1 ] + x[k+  1 *14];
            g[2 ] = u2 * d[2 ] - g[2 ] + x[k+  2 *14];
            g[3 ] = u2 * d[3 ] - g[3 ] + x[k+  3 *14];
        }
        rr[0+14*i] = u * g[0 ] - d[0 ] + x[0 *14] * 0.5;
        rr[1+14*i] = u * g[1 ] - d[1 ] + x[1 *14] * 0.5;
        rr[2+14*i] = u * g[2 ] - d[2 ] + x[2 *14] * 0.5;
        rr[3+14*i] = u * g[3 ] - d[3 ] + x[3 *14] * 0.5;

        d[4 ] = 0; g[4 ] = x[13+14*4 ];
        d[5 ] = 0; g[5 ] = x[13+14*5 ];
        d[6 ] = 0; g[6 ] = x[13+14*6 ];
        d[7 ] = 0; g[7 ] = x[13+14*7 ];
        for (k = 11; k >= 1; k-=2) {
            d[4 ] = u2 * g[4 ] - d[4 ] + x[k+1+4 *14];
            d[5 ] = u2 * g[5 ] - d[5 ] + x[k+1+5 *14];
            d[6 ] = u2 * g[6 ] - d[6 ] + x[k+1+6 *14];
            d[7 ] = u2 * g[7 ] - d[7 ] + x[k+1+7 *14];
            g[4 ] = u2 * d[4 ] - g[4 ] + x[k+  4 *14];
            g[5 ] = u2 * d[5 ] - g[5 ] + x[k+  5 *14];
            g[6 ] = u2 * d[6 ] - g[6 ] + x[k+  6 *14];
            g[7 ] = u2 * d[7 ] - g[7 ] + x[k+  7 *14];
        }
        rr[4+14*i] = u * g[4 ] - d[4 ] + x[4 *14] * 0.5;
        rr[5+14*i] = u * g[5 ] - d[5 ] + x[5 *14] * 0.5;
        rr[6+14*i] = u * g[6 ] - d[6 ] + x[6 *14] * 0.5;
        rr[7+14*i] = u * g[7 ] - d[7 ] + x[7 *14] * 0.5;

        d[8 ] = 0; g[8 ] = x[13+14*8 ];
        d[9 ] = 0; g[9 ] = x[13+14*9 ];
        d[10] = 0; g[10] = x[13+14*10];
        d[11] = 0; g[11] = x[13+14*11];
        d[12] = 0; g[12] = x[13+14*12];
        d[13] = 0; g[13] = x[13+14*13];
        for (k = 11; k >= 1; k-=2) {
            d[8 ] = u2 * g[8 ] - d[8 ] + x[k+1+8 *14];
            d[9 ] = u2 * g[9 ] - d[9 ] + x[k+1+9 *14];
            d[10] = u2 * g[10] - d[10] + x[k+1+10*14];
            d[11] = u2 * g[11] - d[11] + x[k+1+11*14];
            d[12] = u2 * g[12] - d[12] + x[k+1+12*14];
            d[13] = u2 * g[13] - d[13] + x[k+1+13*14];
            g[8 ] = u2 * d[8 ] - g[8 ] + x[k+  8 *14];
            g[9 ] = u2 * d[9 ] - g[9 ] + x[k+  9 *14];
            g[10] = u2 * d[10] - g[10] + x[k+  10*14];
            g[11] = u2 * d[11] - g[11] + x[k+  11*14];
            g[12] = u2 * d[12] - g[12] + x[k+  12*14];
            g[13] = u2 * d[13] - g[13] + x[k+  13*14];
        }
        rr[8 +14*i] = u * g[8 ] - d[8 ] + x[8 *14] * 0.5;
        rr[9 +14*i] = u * g[9 ] - d[9 ] + x[9 *14] * 0.5;
        rr[10+14*i] = u * g[10] - d[10] + x[10*14] * 0.5;
        rr[11+14*i] = u * g[11] - d[11] + x[11*14] * 0.5;
        rr[12+14*i] = u * g[12] - d[12] + x[12*14] * 0.5;
        rr[13+14*i] = u * g[13] - d[13] + x[13*14] * 0.5;

        x += 196;
    }
}

static void _clenshaw_d1(double *rr, const double *x, double u, FINT nroot)
{
    FINT i;
    double d0, d1, g0, g1;
    double u2 = u * 2.;

    for (i = 0; i < nroot-1; i+=2) {
        d0 = 0;
        d1 = 0;
        g0 = x[13+   14*i];
        g1 = x[13+14+14*i];
        d0 = u2 * g0 - d0 + x[12+   14*i];
        d1 = u2 * g1 - d1 + x[12+14+14*i];
        g0 = u2 * d0 - g0 + x[11+   14*i];
        g1 = u2 * d1 - g1 + x[11+14+14*i];
        d0 = u2 * g0 - d0 + x[10+   14*i];
        d1 = u2 * g1 - d1 + x[10+14+14*i];
        g0 = u2 * d0 - g0 + x[9 +   14*i];
        g1 = u2 * d1 - g1 + x[9 +14+14*i];
        d0 = u2 * g0 - d0 + x[8 +   14*i];
        d1 = u2 * g1 - d1 + x[8 +14+14*i];
        g0 = u2 * d0 - g0 + x[7 +   14*i];
        g1 = u2 * d1 - g1 + x[7 +14+14*i];
        d0 = u2 * g0 - d0 + x[6 +   14*i];
        d1 = u2 * g1 - d1 + x[6 +14+14*i];
        g0 = u2 * d0 - g0 + x[5 +   14*i];
        g1 = u2 * d1 - g1 + x[5 +14+14*i];
        d0 = u2 * g0 - d0 + x[4 +   14*i];
        d1 = u2 * g1 - d1 + x[4 +14+14*i];
        g0 = u2 * d0 - g0 + x[3 +   14*i];
        g1 = u2 * d1 - g1 + x[3 +14+14*i];
        d0 = u2 * g0 - d0 + x[2 +   14*i];
        d1 = u2 * g1 - d1 + x[2 +14+14*i];
        g0 = u2 * d0 - g0 + x[1 +   14*i];
        g1 = u2 * d1 - g1 + x[1 +14+14*i];
        rr[i]   = u * g0 - d0 + x[14*i   ] * 0.5;
        rr[i+1] = u * g1 - d1 + x[14*i+14] * 0.5;
    }
    if (i < nroot) {
        d0 = 0;
        g0 = x[13+14*i];
        d0 = u2 * g0 - d0 + x[12+14*i];
        g0 = u2 * d0 - g0 + x[11+14*i];
        d0 = u2 * g0 - d0 + x[10+14*i];
        g0 = u2 * d0 - g0 + x[9 +14*i];
        d0 = u2 * g0 - d0 + x[8 +14*i];
        g0 = u2 * d0 - g0 + x[7 +14*i];
        d0 = u2 * g0 - d0 + x[6 +14*i];
        g0 = u2 * d0 - g0 + x[5 +14*i];
        d0 = u2 * g0 - d0 + x[4 +14*i];
        g0 = u2 * d0 - g0 + x[3 +14*i];
        d0 = u2 * g0 - d0 + x[2 +14*i];
        g0 = u2 * d0 - g0 + x[1 +14*i];
        rr[i] = u * g0 - d0 + x[14*i] * 0.5;
    }
}

static void _matmul_14_14(double *imc, double *im, FINT nroot)
{
    double o7 = 0.14285714285714285714;
    double d0[14];
    FINT i, j;
    for (i = 0; i < nroot; i++) {
        d0[0 ] = 0;
        d0[1 ] = 0;
        d0[2 ] = 0;
        d0[3 ] = 0;
        d0[4 ] = 0;
        d0[5 ] = 0;
        d0[6 ] = 0;
        d0[7 ] = 0;
        d0[8 ] = 0;
        d0[9 ] = 0;
        d0[10] = 0;
        d0[11] = 0;
        d0[12] = 0;
        d0[13] = 0;
        for (j = 0; j < 14; j++) {
            d0[0 ] += im[j+14*i] * COS_14_14[j*14+0 ];
            d0[1 ] += im[j+14*i] * COS_14_14[j*14+1 ];
            d0[2 ] += im[j+14*i] * COS_14_14[j*14+2 ];
            d0[3 ] += im[j+14*i] * COS_14_14[j*14+3 ];
            d0[4 ] += im[j+14*i] * COS_14_14[j*14+4 ];
            d0[5 ] += im[j+14*i] * COS_14_14[j*14+5 ];
            d0[6 ] += im[j+14*i] * COS_14_14[j*14+6 ];
            d0[7 ] += im[j+14*i] * COS_14_14[j*14+7 ];
            d0[8 ] += im[j+14*i] * COS_14_14[j*14+8 ];
            d0[9 ] += im[j+14*i] * COS_14_14[j*14+9 ];
            d0[10] += im[j+14*i] * COS_14_14[j*14+10];
            d0[11] += im[j+14*i] * COS_14_14[j*14+11];
            d0[12] += im[j+14*i] * COS_14_14[j*14+12];
            d0[13] += im[j+14*i] * COS_14_14[j*14+13];
        }
        imc[0 +14*i] = o7 * d0[0 ];
        imc[1 +14*i] = o7 * d0[1 ];
        imc[2 +14*i] = o7 * d0[2 ];
        imc[3 +14*i] = o7 * d0[3 ];
        imc[4 +14*i] = o7 * d0[4 ];
        imc[5 +14*i] = o7 * d0[5 ];
        imc[6 +14*i] = o7 * d0[6 ];
        imc[7 +14*i] = o7 * d0[7 ];
        imc[8 +14*i] = o7 * d0[8 ];
        imc[9 +14*i] = o7 * d0[9 ];
        imc[10+14*i] = o7 * d0[10];
        imc[11+14*i] = o7 * d0[11];
        imc[12+14*i] = o7 * d0[12];
        imc[13+14*i] = o7 * d0[13];
    }
}

void CINTstg_roots(FINT nroots, double ta, double ua, double* rr, double* ww)
{
  const double* x = DATA_X + (nroots-1)*nroots/2 * 19600;
  const double* w = DATA_W + (nroots-1)*nroots/2 * 19600;
  double u, uu, t, tt;
  FINT k, iu, it;
  FINT offset;
  double im [14*nroots];
  double imc[14*nroots];

  t = ta;
  if (t > 19682.99) t = 19682.99;
  u = ua;
  if (t > 1.0) {
      tt = log(t) * 0.9102392266268373 + 1.0; // log(3)+1 
  } else {
      tt = sqrt(t);
  }
  uu = log10(u);

  it = (FINT)tt;
  tt = tt - it;
  tt = 2.0 * tt - 1.0;

  iu = (uu + 7); // 0 <= iu <= 9
  if (iu < 0 || iu > 10) {
      fprintf(stderr, "current implementation assumes 1.0e-7 < U < 1.0e3");
      exit(1);
  } else {
      uu = uu - (iu - 7);
      uu = 2.0 * uu - 1.0;
  }

  offset = nroots * 196 * (iu + it * 10);
  _clenshaw_dc(im, x+offset, uu, nroots);
  _matmul_14_14(imc, im, nroots);
  _clenshaw_d1(rr, imc, tt, nroots);
  _clenshaw_dc(im, w+offset, uu, nroots);
  _matmul_14_14(imc, im, nroots);
  _clenshaw_d1(ww, imc, tt, nroots);
  uu = 1./sqrt(ua);
  for (k = 0; k < nroots; k++) {
      ww[k] *= uu;
  }
}

